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Abstract 

In this paper we present two new methods based on an imphcit Runge-Kutta 
method Gauss which is of algebraic order fourth and has two stages: the first 
one has zero dispersion and the second one has zero dispersion and zero 
dissipation. The efficiency of these methods is measured while integrating 
the radial Schrodinger equation and other well known initial value problems. 
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1. Introduction 

We consider the radial Schrodinger equation: 
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where is the centrifugal potential, V{x) is the potential, E is the 

energy and W{x) = ^^^^ + V{x) is the effective potential. It is valid that 

lim V{x) = 

«— »oo 

and therefore 

lim W{x) = 0. 

We will study the case of E > 0. 

If we divide [0, oo] into small subintervals [oj, bi] so that W{x) is consid- 
ered constant with value Wi, then the problem (1) is reduced to the approx- 
imation 

y'i = (W — E) Ui, whose solution is 

yi{x) = Ai exp (^Vw - Ex^ + Bi exp (^-Vw - Ex^ , Ai, e 3?. (2) 

This form of Schrodinger equation shows why phase fittin is so important 
when new methods are constructed. In the next section we will present the 
most important parts of the theory used. 

The structure of the paper is as follows. Firstly in section 2, the basic 
theory of implicit Runge-Kutta methods is presented. In section 3, the con- 
struction of the methods is introduced. In section 4 and 5, the calculation of 
the algebraic order and the symplecticity respectively of the new methods is 
given. Then in section 6, the numerical results of the new methods are pre- 
sented compared to classical RK methods from the literature while integrat- 
ing well known initial value problems and the radial Schrodinger equation. 
Finally in section 7 our conclusions are presented. 

2. Basic Theory 

2.1. Implicit Method. 

The general form of an s-stage implicit Runge-Kutta method used for the 
computation of the approximate value of yn+i{x) in Problem ([1]), when yn{x) 
is known, is given from the following procedure: 

s 

Wi = f(tn + Cih, yn + hJ2 aijWj) 

(3) 

yn+i = yn + hYj hWi 

i=l 



2 



when at least one a^j 7^ exists with i < j. 

An imphcit Runge-Kutta method can also be presented using the Butcher 
table below: 
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2.2. Phase-Lag Analysis 

Let A and B, s x s matrices, be defined by A 



[1 < i,j < s) and 



bj), (1 <i,j < s) respectively. When the method ([3]) is applied 



B = {a 
to the linear equation 

y' = qy, q eC 
the numerical solution is given by 

det{I - zB) 



yn+i = P{,z)yn, P{z) 



det{I - zA) 



(5) 



(6) 



and can be written in the form 



P{z) = K{v) + iL{v), z = hq (7) 
where K{v) and L{v) are functions of v and i = 

Definition 1. [6] In the implicit s-stage Runge-Kutta method, presented 
in (jlj), the quantities 

(f){v) =v-arg{P{iv)), a{v) = 1 - \P{iv)\, ve^ (8) 

are respectively called the phase-lag or dispersion error and the dissipative 
error. If (j){v) = 0{v'^~^^) and a{v) = O^v"^^^) then the method is said to be 
of dispersive order q and dissipative order r. 



3 



2.3. Stability 

Definition 2. 15| The stability function for an implicit Runge-Kutta method 



is the rational function 

where the vector e = (1, 1)-'", and that a method is A-stable if |-R(-2)| < 1, 
whenever Re{z) < 0, where Re{z) is the real part of z. 

3. Construction of the new Runge-Kutta methods 

We consider the implicit Runge-Kutta method of Gauss, which is of al- 
gebraic order fourth and has two stages. The coefficients are shown in Table 
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Below we present the construction of the methods. 



3.1. Construction of the new method with zero phase lag 

We consider all the values of Table [TUl except 62- By evaluating the phase- 
lag of this method, defined in Definition [H and by solving (f){v) = towards 
62, the result is: 

^ _ 1 -6 + 6 + 72 V - tan{v) v"^ + 2A tan ( v) 

2 tan('u) 

+ tan(i;)i;^y3 - 12tan(i;)t;2y3 - 144tan(t;) ^ ' 

-36 1;2 tan [v] - 12 v'^ tan [v) V3-72v 

The Taylor series expansion of 62 is shown below: 



b2, , = - + —V^ + i — —V?>\v^ + 

^taylor ^ 720 V 6720 8640 
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In the last equation we observe that 



hm 62 



1 



'taylor 



2 



namely when the step-length tends to zero the coefficient of the method 
Gauss appears. 

3.2. Construction of the new method with zero phase lag and zero dissipation 
We consider all the values of Table [TOl except two: 62 and 022- An extra 
equation (apart from the equation of the phase lag) must hold, in order to 
achieve zero phase-lag and zero dissipation. The two equations are = 
and a{y) = 0. 

After satisfying the above two equations, by solving towards 62 and 022, 
the result is: 



bo = - ■ — 

6 B 



(12) 



where 



A 



24 sin {v)vV3 + 2 sin (v) v^Vs + 12 v'^ cos [v) 
36 sin {v)v — 3 sin (v) — 144 cos (v) 





and 



1 C 



(13) 



12 ' D 



with 
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C = 41472 + 72 t;^v^ + 864 90 t;^v^+ 1152?/+ 162 i;^ 
+12096 v^Vs - 14688 (cos {v)f - 144 (cos (t;))^ 
+41472 (cos {v)f - 120 sin (v) V cos (v) 
—864 sin (t;) t;To cos (v) — 504 Tq sin (v) ti^(cos {v))'^ cos (v) 
—v^ sin [v ) v/3cos (t;) + 20736 sin {v) V Vs cos (v) 
-1872 sin (t;) tj'^VScos (v) - 288 Tq (cos (t^))^ v^V^ (cos (t"))^ 
+48z;^ro (cos (t;))^ ^3 (cos {v)f + 576i;2To (cos {v)f (cos (t;))^ 
—36 sin (v) cos (i;) + 10368 sin (v) v cos {v) 
+1728 sin [v) cos {v) -\-?>v'^ sin (i;) cos [v] 
-36 v% (cos {v)f + 132 sin (v) v^V^Tq cos (v) 
+i;^ sin {v) VSTq cos (v) - 1728 sin {v) v VSTq cos {v) 
+18 1;^ (cos {v)f - 6912 Tq (cos {v)f 6 v^VS (cos (z;))^ 
-1728 v'^y/^ (cos (v))^ - 360 v^^/^ (cos (v))^ 



D = -144i;^(cos(t;))^ -6v^(cos(v))^ + 576t;2 ^3 + 864^2 

-144 sin {v)v^^/?,cos (v) - 2 v^v^ - 6 + 120 (cos {v)f^ 

-36 cos {v) sin {v) - 1728 v^Vs (cos (v))^ 

—24 sin (v) v^-\/3cos (v) — 1152 -u^ cos (v) sin {v) 

+3456 cos {v) sin [v) f + 96 -n/S (cos (f))^ — sin (v) -s/S cos (v) 

+24To sin {v) cos {v) + 12 1;^ (cos {v)f 

+ 14 t;^A/3 (cos ("u))^ + 'iv'^ sin (z;) cos (i;) 

sin (t;) cos {v) — 12 sin (i;) v cos (t;) 
-288 sin {v) cos (z;) - 24:V^V3 - 864 v (cos (z;))^ 

where 

^ _ / 2z;4V3-4z;4 + 24?-2-24z;2V3-144 
y (cos(i;))^ 

The Taylor series expansion of 62 and 022 are shown below: 
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1 1 4 1 5^3-8 6 

\ , = - H w H ^ 

^taylor ^ 720 10080 -3 + ^3 



_1 1 5^3-9 4 1 220 73-381 g 
«22,„,,„. - ^ + ^^60 -2 + ^3 ^ ~ 181440 (_2 + 73)^ ^ ^ "' 

In the last equations we observe that the hmits when t; — are equal to 
the corresponding coefficients of the Gauss method. 

4. Algebraic order of the new methods 

The following 8 equations must be satisfied so that the new method main- 
tains the fourth algebraic order of the corresponding classical method pre- 
sented in Table [10 The number of stages is symbolized by s, where s = 2. 



1st Alg. Order (1 equation) 

th = l 

i=l 

2st Alg. Order (2 equations) 

s ^ 
i=l 

3st Alg. Order (4 equations) 

t = I 

i=l 
s 

bittijCj = g 

4st Alg. Order (8 equations) 

t b^cf = \ 

i=l 



(14) 



(15) 



(16) 



bi Ci Ojij Cj g 



i,j=i (17) 

-2 = J_ 
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s ^ 
biClijCljkCk 24 
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4-1. Remainders for the first method (algebraic conditions) 

We present the remainders of the eight equations, that is the difference 
of the right part minus the left part, for the first method: 
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0. This 



We see that the eight equations are held, when h 
means that the new method maintains the algebraic order of the correspond- 
ing classical method. 



4-2. Remainders for the second method (algebraic conditions) 

Now we present the remainders of the equations for the second method: 
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(19) 

We see that for v = the eight equations are held for this method too. Thus 
the new method has also fourth algebraic order. 
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5. Symplecticity of the new methods 

Theorem The Runge-Kutta method (3)- (4) is symplectic when the fol- 
lowing equalities are satisfied 

bittij + bjUji = bibj, I <i,j < s. (20) 

As a classical example we mention the Gauss methods as symplectic 
Runge-Kutta methods. It should be noted that symplectic Runge-Kutta 
methods are always implicit. 

Thus according to the above theorem, the three equations must be satis- 
fied so that the symplecticity of the new methods will be maintained. 

6iaii + biau = bibi (21) 

62^22 + ^2022 = ^2^2 (22) 

6iai2 + 62021 = bib2 (23) 

5.1. Remainders for the first method (simplecticity conditions) 

We present the remainder of the three equations, that is the difference of 
the right part minus the left part, for the first method: 

renii = 

rem, = -j^v'+{-j^ + j^^V^v' + ..^. ^ (24) 

rems = (-3^ + 4^ V^) v ^ + (-241920 + 24TI20 ^) + • • • 

We sec that for the three equations arc held, when h —> => v —> 0. That 
means that the new method maintains the symplecticity of the corresponding 
classical method. 

5.2. Remainders for the second method (simplecticity conditions) 

Now we present the remainders of the equations for the second method: 



rerrii = 

rpm — 1 -45+26 „,4 , 1 545^3-944 6 i , 

rem2- (.2+^3) (_3+y3)^^ + ™ (^^Wif^Wif (25) 

rpm — 1 -5+3 V3 ,, 4 1 -54+31 y/S ,. 6 i 

' '^'"'3 2880 -3+V3 120960 -3+^3 ^ ' ' ' 
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We see that for v = the three equations are held for this method too. 
Thus the new method is also symplectic. 

6. Numerical Results 
6.1. The methods 

In order to measure the efficiency of the methods constructed in this paper 
we compare them to some already known methods, presenting the results of 
the best six. 

I. Method G2-PL-D constructed in this paper, where G2-PL-D means 
the method Gauss two-stages, fourth-order with zero phase-lag and zero dis- 
sipation. 

II. Method G2-PL constructed in this paper, where G2-PL means the 
method Gauss two-stage, fourth-order with zero phase-lag. 

III. Method G2: The classical two-stages and fourth-order Gauss 
method (see [l6|). 

IV. Method SDIRK(3,6,3): The Singly Diagonally-Implicit Runge- 
Kutta method of J. M. Franco, I. Gomez, L. Randez, is third-stage, third 
algebraic order, sixth dispersive order and third dissipative order (see (25|). 

V. Method Radau I: The classical third order Radau method (see 



VI. Method Lobatto IIIC: The classical fourth order Lobatto method 
(see 

6.2. The Problems 

6.2.1. Inverse Resonance Problem 

The efficiency of the two new constructed methods will be measured 
through the integration of problem ([1]) with Z = at the interval [0, 15] 
using the well known Woods-Saxon potential 




(26) 




a 



so for large x 



(asymptotic region) the Schrodinger equation ([T]) becomes 
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y"{x) 



+ V{x) - E 



y{x) 



(27) 



The last equation has two hnearly independent solutions kxji{kx) and 
kxni{kx), where ji and ni are the spherical Bessel and Neumann functions. 
When a; oo the solution takes the asymptotic form 



y{x) ~ Akxji{kx) — Bkxni{kx) (28) 
D [sin{kx - 7r//2) + tan{5i)cos{kx - 7r//2)] (29) 

where 5i is called scattering phase shift and it is given by the following 
expression: 

y[xi+i)C[xi) - y{xi)C{xi+i) 

where S{x) = kxji{kx), C{x) = kxni{kx) and Xj < Xj+i and both belong 
to the asymptotic region. Given the energy we approximate the phase shift, 
the accurate value of which is 7r/2 for the above problem. We will use two 
values for the energy: 989.701916 and 341.495874. As for the frequency w 
we will use the suggestion of Ixaru and Rizea in ^ and 0.: 

r VE^, tf XE[0, 6.5] 

\ ^/E, else x e [6.5, 15] ^ ' 

In Figure 1 we use E = 989.701916 and in Figure 2 we use E = 341.495874. 

6.2.2. Inhomogeneous Equation 

y" = -100y + 99sin{t), with y{0)=l, y'{0)=ll, t e [0, lOOOvr]. Theoretical 
solution: y{x) = sin(t) + sin(lOt) + cos(lOt). 
Estimated frequency: w=10. 

6.2.3. Duffing Equation 

y" = -y - y3 ^ 0.002 cos(l.Olt) , with ?/(0)=0.200426728067, y'{0)=0, 
t e [0,10007r]. Theoretical solution: y{x) = 0.200179477536 cos(l.OU) + 
2.46946143 ^lO"^ cos(3.03t)+3.04014*10-^ cos(5.05t)+3.74*10-i° cos(7.07t) + 

Estimated frequency: w=l. 
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6.2.4- Nonlinear Problem 

y" = -lOOy + sin{y), with y{0)=0, y'{0)=l, t e [0,207r].The theoretical 
solution is not know, but we use ?/(207r) = 3.92823991 * 10~^ and w=10 as 
frequency of this problem. 

7. Conclusions 

In the following figures we present the accuracy of the tested methods 
in connection with loglO of the total stepsx stages of the methods. Firstly, 
in figured] (Resonanse Problem) we use E=989. 701916 and we observe that 
the second method, which has zero phase-lag and zero dissipation, has such 
the same accuracy with the first method, which has zero phased-lag, but 
much better accuracy than the other methods (G2, SD1RK(3,6,3) , Radau I 
, Lobatto IIIC), while in figure [2] (Resonanse Problem) we use E=341. 495874 
and we also see that the second method, which has zero phase-lag and zero 
dissipation, has the same accuracy with the first method, which has zero 
phase-lag, but better accuracy than the other methods (G2, SDIRK(3,6,3) , 
Radau I , Lobatto IIIC). The conclusion from the above is that the difference 
in efficiency is higher when using higher energy. Secondly in figure [3] (Inho- 
mogeneous Equation) it can be observed that the second method developed 
here with zero phase-lag and zero dissipation has way higher efficiency than 
the first method developed here with zero phase-lag, which has far higher ef- 
ficiency than the other methods G2, SDIRK(3,6,3) , Radau I , Lobatto IIIC. 
Thirdly, in figure H] (Duffing Equation) the two new methods have almost 
the same efficiency but much better than the other four methods. Finally in 
figure [5] (Nonlinear Equation) the two new methods have almost the same 
efficiency but much better than the other four methods. This shows the 
importance of zero phase-lag and zero dissipation in this type of problems. 

The stability region for the three methods for v = u h = 50 are shown 
in Figure [61 The stability regions of the methods include the exterior of 
the curves. Note that the curves of all the three methods include the left 
half-plane. 
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ures of the Numerical Results 



Woods-Saxon Potential (E=989.701916) 

9| | ' 1 




steps X stages 



Figure 1: Resonance Problem using £==989. 701916 
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Woods-Saxon Potential (E=341 .495874) 

9| | ' , 




steps X stages 



Figure 2: Resonance Problem using E=341. 495874 
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